Smooth k~temperature models instead of Sharpe-Schoolfield?
Author
Max Lindmark, Jan Ohlberger, Anna Gårdmark
Published
June 15, 2023
Load libraries
Code
pkgs <-c("tidyverse", "tidylog", "sdmTMB", "RColorBrewer", "viridis")## minpack.lm needed if using nlsLM()if(length(setdiff(pkgs, rownames(installed.packages()))) >0){install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T) }invisible(lapply(pkgs, library,character.only = T))#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──#> ✔ dplyr 1.1.2 ✔ readr 2.1.4#> ✔ forcats 1.0.0 ✔ stringr 1.5.0#> ✔ ggplot2 3.4.2 ✔ tibble 3.2.1#> ✔ lubridate 1.9.2 ✔ tidyr 1.3.0#> ✔ purrr 1.0.1 #> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──#> ✖ dplyr::filter() masks stats::filter()#> ✖ dplyr::lag() masks stats::lag()#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors#> #> Attaching package: 'tidylog'#> #> #> The following objects are masked from 'package:dplyr':#> #> add_count, add_tally, anti_join, count, distinct, distinct_all,#> distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,#> full_join, group_by, group_by_all, group_by_at, group_by_if,#> inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,#> relocate, rename, rename_all, rename_at, rename_if, rename_with,#> right_join, sample_frac, sample_n, select, select_all, select_at,#> select_if, semi_join, slice, slice_head, slice_max, slice_min,#> slice_sample, slice_tail, summarise, summarise_all, summarise_at,#> summarise_if, summarize, summarize_all, summarize_at, summarize_if,#> tally, top_frac, top_n, transmute, transmute_all, transmute_at,#> transmute_if, ungroup#> #> #> The following objects are masked from 'package:tidyr':#> #> drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,#> spread, uncount#> #> #> The following object is masked from 'package:stats':#> #> filter#> #> #> Loading required package: viridisLite# devtools::install_github("seananderson/ggsidekick") ## not on CRAN library(ggsidekick)theme_set(theme_sleek())# Check the temperature model script! This is the order based on mean temperature, which makes sense for the sharpe schoolfield plot, and therefore we might keep it across plotsorder <-data.frame(area =c("SI_HA", "BT", "TH", "SI_EK", "FM", "JM", "MU", "FB", "VN", "HO", "BS", "RA")) nareas <-length(unique(order))colors <-colorRampPalette(brewer.pal(name ="RdYlBu", n =10))(nareas)# Load functionshome <- here::here()# fxn <- list.files(paste0(home, "/R/functions"))# invisible(sapply(FUN = source, paste0(home, "/R/functions/", fxn)))
Read and trim data
Code
VBG <-read_csv(paste0(home, "/output/vbg.csv"))#> Rows: 382 Columns: 10#> ── Column specification ────────────────────────────────────────────────────────#> Delimiter: ","#> chr (2): area, type#> dbl (8): cohort, k_mean, k_median, linf_median, k_lwr, k_upr, temp, n#> #> ℹ Use `spec()` to retrieve the full column specification for this data.#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Fit smooth models of k vs temperature
Code
preds <-list()for(a inunique(VBG$area)) { d <-filter(VBG, area == a) m <-sdmTMB(k_mean ~s(temp, k =3), data = d,spatial ="off",spatiotemporal ="off"#,#family = student(df = 6) ) res <-residuals(m)qqnorm(res); qqline(res) nd <-data.frame(area = a,temp =seq(min(d$temp), max(d$temp), length.out =50)) nd$pred <-predict(m, newdata = nd)$estprint(ggplot(nd, aes(temp, pred)) +geom_point(data = d, aes(temp, k_median), inherit.aes =FALSE) +geom_line()) preds <-bind_rows(preds, nd)}#> filter: removed 319 rows (84%), 63 rows remaining